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FINAL  REPORT  ON  “GEOSPATIAL  REPRESENTATION,  ANALYSIS  AND 
COMPUTING  USING  BANDLIMITED  FUNCTIONS”,  AFOSR  GRANT 

FA9550-07-1-0135 


CORY  AHRENS,  GREGORY  BEYLKIN  AND  KRISTIAN  SANDBERG 


Summary 

The  original  goal  of  the  proposal  was  to  address  three  problem  areas: 

(1)  Develop  methods  for  representation  and  accurate  interpolation  of  space-limited  portions  of 
data  using  bases  for  bandlimited  functions. 

(2)  Develop  numerical  integrators  for  solving  systems  of  the  Ordinary  Differential  Equations 
(ODEs)  that  are  based  on  quadratures  for  bandlimited  functions,  rather  than  polynomials. 

(3)  Develop  representation  of  functions  on  the  sphere,  bandlimited  functions  on  spherical 
patches  and  associated  multiresolution  bases. 

We  obtained  new  significant  results  in  all  three  areas  which  were  published  in  3  papers  and  1 
preprint.  Let  us  summarize  these  results. 

(1)  We  have  transferred  code  that  generates  a  local  approximation  to  gravity  models  to  Bran¬ 
don  Jones,  a  Ph.D.  student  in  Aerospace  Department  at  CU,  who  made  some  improvement 
and  fixed  several  bugs.  The  results  of  testing  these  models  have  been  submitted  as  the 
paper  B.  Jones,  G.  Born  and  G.  Beylkin,  “Comparisons  of  The  Cubed  Sphere  Gravity 
Model  with  the  Spherical  Harmonics”  to  Journal  of  Guidance,  Control,  and  Dynamics,  to 
appear  in  Volume  33,  Number  2  (see  http:/ /dx.doi.org/10. 2514/1. 45336).  During  the 
summer  of  2009  Brandon  Jones  implemented  these  localized  gravity  models  for  NASA  at 
Johnson  Space  center. 

We  also  developed  a  new  approach  to  approximating  solutions  of  Laplace’s  equation  satis¬ 
fying  boundary  conditions  using  Gaussians.  Our  new  representation  inherits  a  multireso¬ 
lution  structure  from  the  Gaussian  approximation,  leading  to  fast  algorithms  for  the  eval¬ 
uation  of  the  solutions.  In  the  case  of  the  sphere,  our  approach  provides  a  foundation 
for  a  new  multiresolution  approach  to  evaluating  and  estimating  models  of  gravitational 
potentials  used  for  satellite  orbit  computations.  These  results  have  been  published  as  a 
part  of  the  paper  G.  Beylkin  and  L.  Monzon,  “Approximation  by  exponential  sums  re¬ 
visited”,  Applied  and  Computational  Harmonic  Analysis,  v.28,  pp.  131-149,  2010,  (see 
http:  / /dx.doi.org/10. 1016/j.acha.  2009. 08.011). 

These  results  served  as  a  starting  point  for  an  STTR  grant  involving  two  companies,  Nu- 
merica  and  Omitron.  Also,  local  gravity  models  generated  interest  at  NASA  as  an  approach 
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for  modeling  gravity  of  irregular  shaped  bodies,  e.g.  asteroids.  The  work  in  this  direction 
has  already  started. 

(2)  Current  methods  for  solving  ODEs,  be  that  multistep  or  Runge-Kutta,  are  based  on  polyno¬ 
mial  approximations  of  functions.  However,  both  recent  and  classical  results  indicate  that 
in  many  situations  bandlimited  functions,  rather  than  polynomials,  provide  a  near  optimal 
tool  for  numerical  integration  and  interpolation  of  functions.  Given  recently  developed  tools 
for  computing  with  bandlimited  functions,  our  goal  has  been  to  demonstrate  that,  by  using 
bandlimited  approximations,  we  gain  advantages  in  numerical  solution  of  the  initial  value 
problem  for  the  ordinary  differential  equations  (ODEs).  Specifically,  we  considered  ODEs 
for  orbit  determination  as  a  practical  application  of  our  approach.  Using  the  Gaussian-type 
quadratures  for  bandlimited  functions,  we  developed  ODE  solvers  that  mimic  the  standard 
implicit  Runge-Kutta  methods  with  Gauss-Legendre  nodes.  Similar  to  such  methods,  our 
method  is  A-stable.  Moreover,  we  show  that  in  spite  of  the  approximate  nature  of  our 
quadratures  (generated  for  a  finite  but  arbitrary  precision),  the  integrators  may  easily  be 
made  symplectic.  The  nodes  of  the  quadratures  do  not  concentrate  excessively  near  the 
end  points  thus  allowing  us  to  compute  large  portions  of  an  orbit  at  once.  This,  in  turn, 
allows  us  to  use  initially  a  simplified  gravity  model  (e.g.,  with  only  3  terms)  to  approximate 
a  large  portion  (e.g.,  1/2  of  a  period)  of  an  orbit  by  rapidly  solving  a  system  of  nonlinear 
equations.  We  then  access  the  full  gravity  model  and  evaluate  the  gravitational  force  at 
the  nodes  that,  as  a  result  of  the  initial  low  order  approximation,  are  located  fairly  close  to 
their  correct  positions.  We  then  adjust  the  orbit  without  accessing  the  full  gravity  model. 
This  results  in  an  essentially  correct  trajectory.  At  this  point  we  may  (and  currently  do) 
access  the  full  gravity  model  one  more  time  to  evaluate  the  gravitational  force  and  perform 
another  iteration.  Thus,  we  access  the  full  gravity  model  at  most  twice  per  node  while 
using  a  number  of  nodes  that  is  substantially  lower  than  in  traditional  methods. 

We  have  transferred  a  copy  of  the  code  to  Jack  M.  Van  Wieren,  USAF  AFSPC  for  testing 
and  comparison  with  existing  solvers.  A  draft  paper  is  attached  to  this  report. 

(3)  We  have  developed  a  new  numerical  approach  for  obtaining  quadratures  on  the  sphere  that 
are  invariant  under  the  icosahedral  group.  These  nearly  optimal  quadratures  integrate  all 
(N  +  l)2  linearly  independent  functions  in  a  rotationally  invariant  subspace  of  maximal 
order  and  degree  N.  The  nodes  of  these  quadratures  are  nearly  uniformly  distributed  and 
the  number  of  nodes  is  only  marginally  more  than  the  optimal  ( N  +  l)2/3  nodes.  Using 
these  quadratures,  we  discretize  the  reproducing  kernel  on  a  rotationally  invariant  subspace 
to  construct  an  analogue  of  Lagrange  interpolation  on  the  sphere.  This  representation 
uses  function  values  at  the  quadrature  nodes.  In  addition,  the  representation  yields  an 
expansion  that  uses  a  single  function  centered  and  mostly  concentrated  at  nodes  of  the 
quadrature,  thus  providing  a  much  better  localization  than  spherical  harmonic  expansions. 
We  show  that  this  representation  may  be  localized  even  further.  We  also  describe  two 
algorithms  of  complexity  0(N3)  for  using  these  grids  and  representations.  Finally,  we  note 
that  our  approach  is  also  applicable  to  other  discrete  rotation  groups.  These  results  have 
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been  published,  see  C.  Ahrens  and  G.  Beylkin,  “Rotationally  Invariant  Quadratures  for 
the  Sphere”,  Proceeding  of  Royal  Society  of  London  A,  vol.  465,  pp. 3103-3125,  2009  (see 

http:  / /dx.doi.org/10.1098/rspa.  2009. 0104). 

These  quadratures  form  a  foundation  for  the  development  of  new  multiresolution  grids  on 
the  sphere  with  the  goal  of  replacing  spherical  harmonics  by  localized  representations  that 
are  more  appropriate  for  gravity  modeling  and  other  tasks  where  physics  requires  local 
approximations.  We  have  started  work  on  the  further  localization  of  representations  of 
functions  on  the  sphere  and  developing  multiresolution  constructions  suitable  for  estimation 
of  functions  on  the  sphere. 

Below  we  provide  background  and  motivations  for  these  developments. 

1.  Local  approximation  of  gravity  models  and  a  new  type  of  ODE  solvers 

Space  surveillance,  navigation  of  aircraft  and  missiles  requires  detailed  gravity  and  elevation 
maps.  The  choice  of  representation  of  such  digital  information  affects  efficiency  and  accuracy  of 
numerical  algorithms  that  use  these  data.  It  is  often  highly  desirable  to  combine  heterogeneous 
data  derived  from  sources  with  different  spatial  resolutions,  accuracies,  and  coverages.  Thus,  there 
is  a  need  to  develop  data  representations  and  algorithms  that  are  flexible  and  efficient  to  satisfy 
all  of  these  requirements. 

Many  of  the  current  mathematical  representations  and  associated  numerical  algorithms  were 
developed  in  the  1950s  and  1960s,  at  a  time  when  computers  were  quite  different  from  what  is 
routinely  available  today.  The  difference  is  not  only  in  the  drastically  improved  speed  of  modern 
computers,  but  also  in  the  size  of  available  memory  (RAM),  the  computer  architecture  (e.g., 
parallel  architecture  of  inexpensive  Graphics  Processing  Units  (GPUs)),  and  the  relative  cost  of 
arithmetic  operations.  These  differences,  as  well  as  the  availability  of  recently  developed  numerical 
methods,  require  rethinking  of  the  traditional  approach  to  many  of  these  problems. 

Spherical  harmonic  expansions  have  been  and  still  are  the  preferred  analytical  tool  for  the 
representation  of  gravitational  data.  Among  their  advantages  for  low  resolution  models,  is  the 
relative  simplicity,  rapid  convergence  (although  only  for  smooth  functions),  relatively  low  cost  of 
evaluation  at  each  point,  relatively  low  memory  requirements,  and  the  availability  of  the  Fast 
Fourier  Transform  (FFT)  for  the  acceleration  of  the  azimuthal  part  of  the  associated  transform. 
As  the  gravity  maps  became  more  detailed,  higher  order  spherical  harmonic  expansions  had  to 
be  used  to  represent  them.  Currently,  expansions  of  order  and  degree  360  are  used  in  certain 
situations,  and  those  of  order  and  degree  720  are  being  constructed;  there  is  every  reason  to  expect 
that  the  spatial  resolution  of  the  maps  will  keep  increasing.  As  a  result  of  this  growth,  several 
problems  with  spherical  harmonic  expansions  have  become  obvious,  namely: 

(1)  the  representation  of  data  is  global;  a  local  change  in  the  model  in  any  one  area  requires 
re-calculation  of  the  whole  expansion. 

(2)  evaluating  expansion  of  the  order  and  degree  n  at  a  single  point  costs  0(n2)  operations;  in 
many  situations  this  cost  is  excessive. 
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(3)  typically  the  data  are  needed  within  a  relatively  small  window,  covering  only  a  small  part 
of  the  sphere;  extracting  these  data  from  the  spherical  harmonic  expansion  is  inconvenient 
and  could  be  expensive. 

Several  local  approximation  of  gravity  models  have  been  developed  to  address  some  of  these  issues. 
In  this  project  we  demonstrated  that  these  approximations  do  not  compromise  accuracy  in  orbit 
determination  but  drastically  improve  the  speed  of  these  computations.  On  the  other  hand,  we 
still  do  not  have  a  fully  adequate  “replacement”  for  spherical  harmonics  and  our  effort  has  been 
advancing  with  this  goal  in  mind. 

A  more  detailed  accuracy  comparison  of  using  local  gravity  models  in  orbit  determination  may 
be  found  in  http:/ /dx.doi.org/10. 2514/1. 45336.  A  new  approach  to  estimating  local  gravity 
models  has  been  suggested  in  http:/ /dx.doi.org/10.1016/j.acha.2009.08.011. 

Other  aspects  of  these  type  of  problems,  namely,  a  new  type  of  A-stable,  symplectic  integrator 
is  describes  in  the  attached  preprint. 

We  also  include  below  a  brief  description  and  motivation  for  constructing  quadratures  invariant 
under  the  icosahedral  group  and  refer  to  http://dx.doi.org/10.1098/rspa.2009.0104  for  a 
more  detailed  description. 

2.  Nearly  optimal  quadratures  on  the  sphere 

Many  problems  in  physics,  mathematics  and  engineering  involve  integration  and  interpolation  on 
the  sphere  in  M3.  Of  particular  importance  are  discretizations  of  rotationally  invariant  subspaces 
of  L 2  (§2)  that  integrate  all  spherical  harmonics  up  to  a  fixed  order  and  degree.  A  typical  approach 
to  discretizing  the  sphere  is  that  of  equally  spaced  discretization  in  azimuthal  angle  and  Gauss- 
Legendre  discretization  in  polar  angle,  leading  to  an  unreasonably  dense  concentration  of  nodes 
near  the  poles.  It  is  well  known  that  in  a  variety  of  applications  such  concentration  of  nodes  may 
lead  to  problems  when  using  these  grids. 

We  develop  a  systematic  numerical  approach  for  constructing  nearly  optimal  quadratures  in¬ 
variant  under  the  icosahedral  group  to  integrate  rotationally  invariant  subspaces  of  L2  (§2)  up  to 
a  fixed  order  and  degree.  Using  these  grids  and  a  reproducing  kernel,  we  show  how  to  replace 
the  standard  basis  of  spherical  harmonics  on  a  rotationally  invariant  subspace  by  a  representation 
formed  using  a  single  function  centered  at  the  quadrature  nodes.  The  reproducing  kernel  is  mostly 
concentrated  near  the  corresponding  grid  point.  In  the  resulting  representation,  the  coefficients, 
up  to  a  factor,  are  the  values  on  the  grid  of  the  function  being  represented.  We  may  interpret  this 
construction  as  an  analogue  of  Lagrange  interpolation  on  the  sphere  and  note  that  it  allows  us  to 
develop  well  conditioned  linear  systems  for  interpolation  in  contrast  to  some  earlier  constructions. 

An  alternative  to  spherical  harmonics  has  long  been  sought,  especially  for  numerical  purposes. 
Spherical  harmonics  provide  an  efficient  orthonormal  basis,  nicely  subdivided  into  rotationally 
invariant  subspaces.  However,  the  global  support  of  these  functions  poses  a  serious  difficulty  in 
problems  where  physical  effects  are  localized.  In  fact,  the  global  nature  of  spherical  harmonics  is 
a  consequence  of  their  optimality.  Therefore,  if  we  want  localized  functions  to  represent  the  same 
subspaces,  we  necessarily  must  have  a  less  efficient  representation. 
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FIGURE  2.1.  Positions  of  7212  quadrature  nodes  of  a  quadrature  integrating  exactly 
all  spherical  harmonics  in  the  subspace  of  maximal  order  and  degree  145.  This 
quadrature  has  efficiency  i /  =  0.98521....  (efficiency?/  =  1  implies  that  each  node 
accommodates  3  degrees  of  freedom  in  the  subspace). 


We  view  our  approach  as  the  first  step  in  constructing  a  local  and  multiresolution  representation 
of  functions  on  the  sphere  that  respects  rotationally  invariant  subspaces.  We  note  that  the  high 
efficiency  of  quadratures  constructed  in  this  paper  implies  a  near  uniform  distribution  of  nodes 
on  the  sphere.  On  the  other  hand,  the  nodes  maintain  a  regular  organization  visually  similar  to 
that  of  geodesic  or  equal  area  grids.  Moreover,  our  grids  are  associated  with  rotationally  invariant 
subspaces,  an  important  property  in  a  number  of  numerical  applications,  e.g.,  geodesy.  To  date, 
we  have  constructed  grids  which  integrate  subspaces  of  maximum  order  and  degree  ranging  from 
5  up  to  210.  As  an  example,  we  illustrate  Figure  2.1  a  grid  with  7212  nodes  integrating  subspaces 
with  maximal  order  and  degree  145. 

The  rotationally  invariant  spherical  grids  constructed  here  have  many  applications.  Let  us 
mention  a  few  specific  problems  in  some  detail.  First,  due  to  the  central  role  played  by  spherical 
harmonics  in  the  theory  of  gravity  and  magnetic  fields,  solutions  to  many  geodetic  problems  use 
them  as  a  basis.  Yet,  their  global  support  is  inconsistent  with  the  physical  nature  of  the  problem 
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leading  to  many  difficulties  in,  e.g.,  constructing  gravity  models.  The  grids  developed  in  this 
paper  provide  a  first  step  toward  replacing  spherical  harmonics  with  localized  functions.  We  plan 
to  continue  work  in  this  direction.  Second,  the  equations  used  in  global  atmospheric  modeling 
are  typically  posed  on  the  sphere.  Current  spectral  methods  which  use  spherical  harmonics  suffer 
from  the  above  mentioned  problems  of  nodal  clustering  and  require  additional  steps  to  alleviate  the 
problem.  The  new  representations  developed  in  this  paper  eliminate  clustering  and  singularities 
due  to  the  coordinate  system  and  should  provide  efficient  solution  methods.  Third,  acoustic  and 
electromagnetic  scattering  problems  posed  as  integral  equations  involve  integration  over  spherical 
domains.  New  algorithms  for  the  numerical  solution  of  these  integral  equations  may  be  based  on 
the  results  of  this  paper.  Finally,  we  mention  a  numerical  technique  used  in  molecular  dynamics 
calculations  known  as  discrete  variable  representation  (DVR).  Our  quadratures  should  extend  such 
methods  by  allowing  effectively  an  arbitrary  order  and  degree. 

Current  address:  Department  of  Applied  Mathematics,  University  of  Colorado  at  Boulder,  526  UCB  Boulder, 
CO  80309-0526, 


ODE  Solvers  Using  Bandlimited 
Approximations 


G.  Beylkin  and  K.  Sandberg 

Department  of  Applied  Mathematics 
University  of  Colorado  at  Boulder 
526  UCB 

Boulder,  CO  80309-0526 


Abstract 

We  use  generalized  Gaussian  quadratures  for  exponentials  to  develop  a  new  ODE 
solver.  These  generalized  Gaussian  quadratures  integrate  functions  elbx  for  all  |6|  <  c, 
where  the  nodes  and  weights  are  computed  for  a  given  bandlimit  c  and  any  user 
selected  accuracy  e.  An  imprortant  property  of  nodes  of  these  quadratures  is  that 
they  do  not  concentrate  excessively  near  the  end  points  of  the  interval  as  the  nodes 
of  the  standard  polynomial-based  Gaussian  quadratures.  The  new  ODE  solver  is 
analogous  to  the  usual  implicit  Runge  Kutta  (collocation)  method  but  it  allows  us 
to  use  a  large  number  of  nodes  due  to  properties  of  the  quadrature.  We  show  that  the 
resulting  ODE  solver  is  symplectic  and  A-stable.  We  use  this  solver  in  the  problem  of 
orbit  determination  and  achieve  speed  close  to  that  of  an  explicit  multistep  method. 


1  Introduction 

Current  methods  for  solving  ODEs,  be  that  multistep  or  Runge-Kutta,  are 
based  on  polynomial  approximations  of  functions.  On  the  other  hand,  both 
recent  and  classical  results  [5,4,18,16,11,17]  indicate  that  in  many  situations 
bandlimited  functions  provide  a  near  optimal  tool  for  numerical  integration 
and  interpolation  of  functions.  For  example,  a  time  domain  solver  for  the  wave 
equation  [5]  uses  bandlimited  approximations  and  yields  about  12  digits  of  ac¬ 
curacy  with  only  3  nodes  per  wavelength.  Given  recently  developed  tools  for 
computing  with  bandlimited  functions,  our  goal  is  to  demonstrate  that,  by  us¬ 
ing  bandlimited  approximations,  we  also  gain  advantages  in  numerical  solution 

1  This  research  was  partially  supported  by  AFOSR  grant  FA9550-07-1-0135,  NSF 
grant  DMS-0612358,  DOE/ORNL  grants  4000038129  and  DE-FG02-03ER25583. 
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of  the  initial  value  problem  for  the  ordinary  differential  equations  (ODEs).  As 
an  example,  we  use  ODEs  for  orbit  determination  to  provide  us  with  both,  a 
practical  application  as  well  as  a  gauge  to  ascertain  the  performance  of  new 
algorithms. 

Numerical  solution  of  ODEs  is  a  mature  area  of  applied  mathematics  with 
many  well-developed  software  packages.  In  spite  of  this  satisfactory  state  of 
affairs,  we  challenge  the  usual  approach  to  selecting  an  ODE  solver  for  a  given 
problem.  We  note  a  specific  inefficiency  of  multistep  methods,  the  fact  that 
such  methods  can  only  be  A-stable  if  their  order  does  not  exceed  2  (the  so- 
called  Dahlquist  barrier).  Whereas  implicit  Runge-Kutta  methods  with  Gauss- 
Legendre  nodes  are  A-stable  and  symplectic,  we  are  limited  in  using  a  large 
number  of  Gauss-Legendre  nodes  per  time  step  since  such  nodes  concentrate 
excessively  near  the  end  points  of  the  time  interval.  We  note  that  unlike  in 
problems  of  wave  propagation  [5]  where  solutions  are  naturally  expressed  via 
exponentials,  solutions  of  some  ODEs  may,  in  fact,  be  polynomials  or  other 
functions  that  do  not  have  an  efficient  approximation  via  exponentials.  For 
such  equations  it  actually  may  be  advantageous  to  use  polynomial  quadra¬ 
tures.  However,  it  is  more  typical  to  encounter  ODEs  where  solutions  do  have 
an  efficient  approximation  via  exponentials  and  our  method  is  most  suitable 
for  such  ODEs. 

Unlike  the  classical  Gaussian  quadratures  for  polynomials,  the  Gaussian  type 
quadratures  for  exponentials  attempt  to  integrate  an  infinite  set  of  functions, 
namely,  elbx  with  \b\  <  c,  using  a  finite  set  of  nodes.  Clearly,  there  is  no  way  to 
accomplish  this  exactly.  Thus,  these  quadratures  are  constructed  so  that  all 
exponentials  for  \b\  <  c  are  integrated  with  accuracy  of  at  least  e,  where  e  is 
arbitrarily  small  but  finite.  Such  quadratures  were  constructed  in  [4,18]  (we  use 
those  in  [4]).  An  important  observation  in  [5]  is  that  the  nodes  of  quadratures 
of  this  type  do  not  concentrate  excessively  near  the  end  points,  thus  allowing 
us  to  use  as  many  nodes  as  necessary,  without  a  penalty  due  to  their  number. 
The  density  of  nodes  increases  toward  the  end  points  of  the  interval  only  by 
a  factor  that  depends  on  the  desired  accuracy  but  not  on  the  overall  number 
of  nodes.  Since  we  may  choose  many  nodes,  say  M,  it  makes  sense  to  ask 
if  the  integration  matrix  of  an  implicit  Runge-Kutta  type  method  can  be 
applied  in  0(M )  or  0(M  log  M  operations  rather  0(M2).  We  note  that  such 
question  does  not  make  sense  for  the  classical,  polynomial-based  quadratures 
since  a  rapid  concentration  of  nodes  towards  the  end  points  severely  limits  the 
number  of  nodes.  We  show  that,  indeed,  it  is  possible  to  apply  the  integration 
matrix  in  a  fast  manner  (within  a  finite  accuracy  e).  This  may,  for  example, 
be  accomplished  by  using  the  partitioned  low  rank  (PLR)  representation  as  it 
was  described  in  [BEY-SAN:2005].  Using  PLR  representation,  the  speed  of  an 
implicit  method  is  now  of  the  same  order  as  that  of  an  explicit  method.  We 
note  that  there  might  be  even  more  efficient  approaches  to  the  fast  application 
of  the  integration  matrix  than  PLR  representation. 
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This  change  of  complexity  allows  us  to  challenge  the  usual  view  that  explicit 
methods  are  inherently  faster  than  implicit  methods.  While  explicit  methods 
require  a  significant  oversampling,  implicit  methods  only  needs  to  maintain 
minimally  appropriate  sampling  but  do  extra  work  to  solve  a  system  of  equa¬ 
tions  at  each  step.  If  the  efficiency  is  measured  in  terms  of  the  number  of 
function  calls,  then  their  total  number  is  no  longer  differs  substantially  be¬ 
tween  explicit  and  implicit  methods. 

Using  the  Gaussian- type  quadratures  for  bandlimited  functions,  we  construct 
ODE  solvers  that  mimic  the  standard  implicit  Runge-Kutta  methods  with 
Gauss- Legendre  nodes  (see  e.g.  |ISERLE:1996]).  Similar  to  such  methods, 
our  method  is  A-stable.  Although  our  algorithms  are  based  on  approximate 
quadratures  and,  thus,  cannot  produce  results  with  accuracy  better  than  the 
chosen  accuracy  e,  we  note  that  if  e  ~  1CD16  the  effect  of  such  limitation  is  the 
same  as  using  double  precision  in  floating  point  arithmetic.  Moreover,  we  show 
that  in  spite  of  the  approximate  nature  of  our  quadratures,  the  integrators  may 
be  made  exactly  symplectic. 

All  these  properties  make  our  approach  attractive  for  a  number  of  applications 
and  we  select  orbit  determination  as  an  example.  Not  only  does  it  provide  us 
with  a  practical  application,  it  also  allows  us  to  demonstrate  certain  “tricks” 
associated  with  our  choice  of  the  method  that  further  reduce  the  computa¬ 
tional  cost.  Although  orbit  determination  typically  is  not  a  stiff  problem,  we 
demonstrate  advantages  of  using  an  implicit  method  in  that  we  substantially 
reduce  the  required  number  of  function  calls  to  the  full  gravity  model  (beyond 
the  reduction  of  overall  number  of  function  calls  in  comparison  with  existing 
methods) . 


2  Preliminaries:  quadratures  for  bandlimited  functions 

2.1  Bandlimited  functions  as  a  replacement  of  polynomials 

Recently  constructed  quadratures  [18,4]  address  efficiency  of  sampling  for  rep¬ 
resenting  functions  by  breaking  with  the  conventional  approach  of  using  poly¬ 
nomials  as  the  fundamental  tool  in  analysis  and  computation.  The  approach 
based  on  polynomial  approximations  has  a  long  tradition  and  leads  to  such 
notions  as  the  order  of  convergence  of  numerical  schemes,  Gaussian  quadra¬ 
tures  (for  polynomials),  polynomial  based  interpolation,  and  so  on.  Recently, 
the  dominance  of  such  an  approach  (it  clearly  remains  reasonable  for  many 
problems)  has  been  successfully  challenged,  ft  turns  out  that  constructing 
quadratures  for  bandlimited  functions,  e.g.,  exponentials  etbx,  with  \b\  <  c, 
where  c  is  the  bandlimit,  a  fixed  parameter,  leads  to  significant  improvement 


3 


in  performance  of  algorithms  for  interpolation,  estimation  and  solving  partial 
differential  equations  [5]. 


2.2  Bases  for  bandlimited  functions 

Whenever  measurements  are  performed,  all  sensors  are  of  limited  size;  their 
frequency  response  for  all  practical  purposes  is  also  limited  outside  a  finite 
range.  On  the  other  hand,  it  is  well-known  that  a  function  whose  Fourier 
Transform  has  compact  support  can  not  have  compact  support  itself  (unless  it 
is  identically  zero).  Therefore,  it  is  natural  to  analyze  an  operator  whose  effect 
on  a  function  is  to  truncate  it  both  in  the  original  and  the  Fourier  domains. 
This  has  been  the  topic  of  a  series  of  seminal  papers  by  Slepian  et  al.  [17], 
[11],  [12],  [14],  [15],  where  it  is  observed  ( inter  alia )  that  the  eigenfunctions 
of  such  operator  (see  (1)  below)  are  the  Prolate  Spheroidal  Wave  Functions 
(PSWFs)  of  classical  Mathematical  Physics. 

While  periodic  bandlimited  functions  may  be  expanded  into  Fourier  series, 
and  band- limited  functions  on  the  real  line  may  be  represented  via  the  Fourier 
Integral  Transform,  we  must  also  deal  with  non-periodic  functions  on  intervals , 
where  neither  the  Fourier  series  nor  the  Fourier  Integral  can  be  used  efficiently. 
This  motivates  the  introduction  of  a  basis  that  efficiently  represents  (on  an 
interval)  functions  of  the  form  elbx  for  an  arbitrary  real  value  b,  as  long  as 
|6|  <  c,  with  a  fixed  (bandlimit)  parameter  c.  Since  b  varies  continuously,  such 
a  basis  is  not  finite.  On  the  other  hand,  an  arbitrary  but  finite  precision  is 
achievable  with  a  finite  basis  consisting,  for  example,  of  appropriately  chosen 
PSWFs. 

For  a  real  number  c  >  0  (to  be  referred  to  as  the  bandlimit),  we  consider  the 
operator  Fc  :  L2[— 1, 1]  — ►  L2[— 1, 1],  defined  by  the  formula 

Fc{i))[u)  =  J  ^  elcxu'i/j(x)dx,  (1) 

and  the  operator  Qc  =  j-F*Fc",  it  is  easily  seen  that 

^  ,  (W  .  1  f1  sin (c(y  —  x))  ,  ,  ,  , 

Qc{iJj){y)  =  -  - — V’O)  dx.  (2) 

7 r  J-i  y  —  x 

The  eigenfunctions  ^q,  ^  •  •  •  of  Qc  coincide  with  those  of  Fc,  and  the 

eigenvalues  fij  of  Qc  are  connected  with  the  eigenvalues  A  j  of  Fc  by  the  formula 

hi  =  )  (3) 

for  all  j  =  0, 1,  2,  •  •  • ,  where  the  eigenvalues  are  ordered  by  decaying  magni¬ 
tude. 
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In  many  respects,  PSWFs  are  strikingly  similar  to  orthogonal  polynomials; 
they  are  orthonormal,  constitute  a  Chebychev  system,  and  admit  a  version 
of  Gaussian  quadratures  (see  [18],  [4]).  One  of  the  results  in  [18]  and  [4]  is 
formulated  as 


Proposition:  For  c  >  0  and  e  >  0,  there  exist  nodes  -1<9i<92<---< 
9m  <  1  and  corresponding  weights  Wk  >  0,  such  that  for  any  x  G  [—1, 1], 


M 


Actx 


' -1 


dt  —  ^2  wke1' 


icOkX 


k= 1 


<  e, 


(4) 


where  the  number  of  nodes,  M,  is  (nearly)  optimal.  The  nodes  and  weights 
maintain  the  natural  symmetry,  9k  =  —  0M-k+i  and  Wk  =  WM-k+i-  When  the 
functions  ip q,  V’i,  ip  2,  ' ' '  >  'P’m- i  are  used  as  a  basis  for  the  interpolation  on  the 
interval  [—1,1],  with  the  points  9\,  92,  ■  •  ■  ,9  m  used  as  the  interpolation  nodes, 
the  resulting  interpolation  formula  is  stable. 

This  proposition  provides  a  tool  for  the  numerical  integration  and  interpo¬ 
lation  of  functions  of  the  form  elbx  on  [—1,1],  The  nodes  and  weights  are 
functions  of  both  the  bandlimit  c  >  0  and  the  accuracy  e  >  0,  and  may  be 
viewed  as  the  generalized  Gaussian  quadratures  for  the  bandlimited  functions. 
It  is  worth  noting  that  the  algorithm  in  [4|  identifies  the  nodes  of  the  gener¬ 
alized  Gaussian  quadratures  in  (4)  as  zeros  of  the  discrete  prolate  spheroidal 
wave  functions  (DPSWFs)  corresponding  to  small  eigenvalues. 

One  of  problems  associated  with  the  numerical  use  of  orthogonal  polynomials, 
is  the  concentration  of  their  roots  near  the  ends  of  the  interval  of  their  defini¬ 
tion.  Using  nodes  9 l7  92,  ■  ■  ■  ,  9m,  we  maintain  a  nearly  optimal  sampling  rate, 
close  to  the  Nyquist  rate  for  periodic  functions.  Let  us  consider  the  ratio 


r(M,  e)  = 


92  -  9\ 


LM/2J 
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[M/2J-1 


(5) 


where  “ \Mj 2j”  denotes  the  least  integer  part.  Observing  that  the  distance  be¬ 
tween  nodes  of  Gaussian  quadratures  changes  monotonically  from  the  middle 
of  an  interval  toward  its  end  points,  and  that  the  smallest  distance  occurs  be¬ 
tween  the  nodes  closest  to  the  end  point,  the  ratio  5  may  be  used  as  a  measure 
of  node  accumulation.  For  example,  the  distance  between  the  nodes  near  the 
end  points  of  the  standard  Gaussian  quadratures  for  polynomials  decreases 
as  0(l/M 2),  so  that  we  have  r(M,  e)  =  (D(l/M),  where  M  is  the  number  of 
nodes, 

In  Figure  1,  considering  bandlimit  c  as  a  function  of  the  number  of  nodes,  M , 
and  the  desired  accuracy  e,  we  observe  that  the  oversampling  factor, 


a(M,  e)  = 


nM 


>  1, 
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Figure  1.  The  ratio  r(M,  e)  in  (5)  and  the  oversampling  factor  a(M,  e)  plotted  against 
the  number  of  nodes  for  quadratures  of  accuracy  e  ~  10^'and  e  ~  lO^1' . 

approaches  1  for  large  M.  This  factor  compares  the  critical  rate  of  sampling 
of  a  smooth  periodic  function,  either  for  integration  or  interpolation,  to  that 
of  smooth  non-periodic  function  defined  on  an  interval.  We  recall  that  in  the 
case  of  the  Gaussian  quadratures  for  polynomials,  the  oversampling  factor 
approaches  ^  rather  than  1  (see  e.g.  [8]). 

To  illustrate  the  gain  in  efficiency,  we  compare  accuracy  of  differentiation  using 
bases  for  bandlimited  functions,  finite  differences  and  spectral  differentiation 
using  the  Chebyshev  polynomials.  We  construct  two  derivative  matrices  us¬ 
ing  bandlimited  bases  with  accuracy  e  =  10“ '  and  bandlimit  c  =  237T,  and 
with  e  =  ICG13  and  bandlimit  c  =  18.57T.  For  comparison,  we  construct  a 
second-order  central  finite-difference  derivative  matrix  and,  for  spectral  differ¬ 
entiation,  a  block  diagonal  spectral  derivative  matrix  with  4  blocks  where  each 
diagonal  block  is  a  derivative  matrix  with  respect  to  the  first  16  Chebyshev 
polynomials.  We  differentiate  the  function  f(x)  =  sin (bx)  for  200  values  of  b , 
— 327T  <  b  <  32tt.  In  Figure  2  we  illustrate  the  extra  bandwidth  where  the 
bandlimited  derivatives  are  accurate  [5]. 


3  Discretization  of  Picard  integral  equation 


Consider  the  initial  value  problem 

y'  =  y(0)  =  y0, 

or,  equivalently, 


y (t)  =  yo  +  /  f(s,y(s))  ds. 
Jo 


(6) 


6 


Figure  2.  Comparison  of  absolute  errors  for  the  first  derivative  of  the  function  sin  (for) 
in  the  interval  [—1,1]  with  |6|  <  327T.  In  all  examples  the  derivative  matrices  use  64 
independent  parameters. 

For  a  given  accuracy  e,  we  select  Gaussian  nodes  {rJ  }jl£1  such  that  on  [0,  £] 

M 

llf(f>y(f))  -  £f(ri>y(ri))^WII  <  c  (7) 

j=l 

where  Rj(t)  are  interpolating  basis  functions  associated  with  the  generalized 
Gaussian  nodes  for  exponentials  [4,5].  We  discretize  (6)  by  using  (7)  to  obtain 
a  nonlinear  system, 


m  T.  M 

y 0;)  =  y0  +  Ef(rJ>y(ri))  /  Rj(s)ds  =  yo  +  '52Sijf(Tj,y(rj)),  (8) 

■ _ i  J  0  ■ _ i 


3= 1  ^  3=1 

where  Stl  =  /0r‘  Rj(s)ds  is  the  integration  matrix. 


3.1  Relation  to  Implicit  Runge  Kutta  methods  based  on  collocation 


Consider  (8)  with  M  quadrature  nodes  (and  the  corresponding  weights 

{wj}1-^)  in  the  interior  of  the  interval  [0,  t\.  Then  implicit  Runge-Kutta  meth¬ 
ods  based  on  collocation  discretize  (6)  as 

M 

y  (t)  =  y0  +  J2  wArv  y(ri))>  (9) 

3= 1 

where  {y(r3)}jll1  are  obtained  by  solving  the  nonlinear  system  (8).  The  nodes, 
weights,  and  the  entries  of  the  integration  matrix  are  typically  organized  in 
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the  Butcher  tableau, 


T 

s 

wl 

Unlike  in  the  standard  implicit  Runge-Kutta  method  based  on  Gauss-Legendre 
quadratures,  we  solve  (8)  on  a  time  interval  containing  a  large  number  of 
quadrature  nodes,  M  (since  these  nodes  do  not  concentrate  excessively  near 
the  end  points). 

3.2  Exact  Linear  Part 

In  many  problems  (including  that  of  orbit  determination),  the  right  hand  side 
of  the  ODE,  f(t,  y),  may  be  split  into  a  linear  and  nonlinear  parts, 

f&y  (*))  =  Ly(t)  +  g(t,y(t))., 

so  that  the  integral  equation  (6)  may  be  written  as 

y(t)=etLy0+[  e(t_s)Lg(s,  y(s))  ds.  (10) 

Jo 

If  the  operator  etL  is  computed  efficiently,  this  formulation  may  lead  to  savings 
when  solving  the  integral  equation  iteratively.  We  discretize  (10)  by  using  (7) 
and  obtain 

M 

y(Ti)  =  eTiLy0  +  ekrTj)Lg(rj,  y (t,-)) 

j= 1 

M 

=  eTiLyo  +  5|e(T]|7i')Lg(rj,  y^-))  (11) 

3= 1 

where  =  /0Tl  Rj(s)ds.  We  note  that  (8)  is  a  special  case  of  (11)  with  L  =  0 
and  g  =  f. 

3.3  Algorithm 

Next  we  describe  a  fixed  point  iteration  to  solve  (11).  Let  Nlt  denote  the 
number  of  iterations,  which  can  either  be  set  to  a  fixed  number  or  deter¬ 
mined  adaptively.  Labeling  the  intermediate  solutions  in  the  iteration  scheme 
as  yW,  i  =  1, ... ,  Nit,  we  have 


Rj(s)ds 


(1)  Initialize  2/(1)(t;)  =  y0,  i  = 

(2)  For  k  =  1, . . . ,  Nit 


For  m  —  1, . . . ,  M 


(a)  Update  the  solution  at  the  node  m: 


(b)  Update  the  right  hand  side:  g (Tm,  (rm)) 


4  Symplectic  integrators 

Following  [13],  let  us  introduce  matrix  M  =  {'nikjjk! j=\  f°r  an  hnplicit  M- stage 
Runge-Kutta  (RK)  scheme, 


(12) 


Wlkj  R kj  T  WjSjk  WpWj , 


where  w  =  {wk]iLi  and  S  =  {Skj}i j=i  dehne  the  Butcher’s  tableau  for  the 
method. 

It  is  shown  in  [13]  that 

Theorem  1  If  matrix  M  =  0  in  (12),  then  an  implicit  M-stage  RK  scheme 
is  symplectic. 

This  condition,  M  =  0,  is  satisfied  for  the  Gauss-Legendre  RK  methods,  see 
e.g.  [6,13].  We  will  enforce  it  to  construct  symplectic  integrators  for  collocation- 
type  methods  based  on  the  generalized  Gaussian  quadratures  for  bandlimited 
functions.  In  spite  of  the  fact  that  such  quadratures  are  accurate  only  up  to 
a  fixed  (but  arbitrary)  accuracy  e,  we  can  still  satisfy  the  condition  M  =  0 
exactly. 

Our  proof  is  based  on  several  observations. 

Proposition  2  Let  {i?i(r)}^1  be  interpolating  basis  functions  for  the  quadra¬ 
ture  nodes  { r, }£L1;  such  that  R-fjj)  =  .  Define  Ffr)  =  ffl  Ri(s)  ds  and  let 

{w]}!]=\  be  quadrature  nodes  such  that 


L 


i 


M 


-1 


CiMC'M  dr  -  J2wiFAT*)Fi(T*)  < j2 


(13) 


k= 1 
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Then 


i'l  1  JZi  Rj(s)  ds  Ri(r)  dr 


i-i 


Rj(s )  ds  — 


Wi 


<  e 


PROOF.  We  first  observe  that  due  to  the  interpolating  property  of  Ri(r), 
we  have 


J!i^(s)  ds  WiRiidi)  _  Ek=i  J-lRjis)  ds  wkRi(9k ) 


Rifs')  ds  = 


i-i 


Wi 


Wi 


(14) 


Next  we  note  that 

pi  rr 


J  J  Rj(s )  ds  Ri(r)  dr  =  j  R,(r)F/(r)  dr, 
so  that  proposition  follows  by  combining  (14)  and  (13). 


Proposition  3  Let  {Rj(r)}f; f1  be  interpolating  basis  functions  for  the  quadra¬ 
ture  nodes  such  that  Ri{jj)  =  5ij.  Let  {wj}jL1  be  the  corresponding 

quadrature  weights  so  that 


i  M 

/  Ri(r)  dr  -  V  wkRi{rk ) 

J~l  k=i 


<  e2. 


Then 


J  Rifs)  ds  —  Wi 


<  e 


PROOF.  The  result  follows  from  the  interpolating  property  of  Ri(r). 

Theorem  4  Given  M  quadrature  nodes  {rk}^l  and  interpolating  functions 
{Rk}kLi,  let  weights  for  the  quadrature  be  defined  as 


<l.‘k  -  j  Rk{r)dr 


(15) 


and  the  integration  matrix  as 


Sm  =  &(&*(')*>)  KM*  M  = 


Wk 


Then 


wkSki  +  wiSik  -  wkwi  =  0, 

and  the  implicit  scheme  using  these  nodes  and  weights  is  symplectic. 
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PROOF.  Using  propositions  above,  we  observe  that  the  weights  defined  in 
(15)  are  the  same  (up  to  accuracy  e2)  as  those  of  the  quadrature.  The  result 
follows  if  we  set  Fk{r)  =  JZ]  Rk(r)  dr,  F[(t)  =  Rk(j)  and  integrate  by  parts 
to  obtain 

wkSki  +  wiSik  -  wkwi  =  J  ^  Fi(r)F'k{r)  dr  +  J  ^  Ffc(r)F/(r)  dr  -  wkwi 

=  Fi(l)Fk(l)  -wkwh 

By  the  definition  of  the  weights  we  have  Fk(  1)  =  Wk  and,  hence,  Fi(l)Fk(l)  — 

wkwi  =  0. 

Remark  5  Implicit  Runge-Kutta  methods  based  on  Gauss- Legendre  nodes  are 
A-stable  (see  e.g  [10]).  We  verified  numerically  (within  the  accuracy  of  under¬ 
lying  quadrature)  that  this  property  carries  on  to  the  schemes  based  on  the 
generalized  Gaussian  quadratures  for  exponentials.  Currently,  we  do  not  have 
a  proof  of  this  observation. 


5  Fast  application  of  integration  matrix 

The  repeated  application  of  the  integration  matrix  S  dominates  the  com¬ 
putational  cost  of  current  algorithm  since  applying  S'  as  a  dense  matrix  re¬ 
quires  0(M2)  operations.  We  demonstrate  that  this  cost  may  be  reduced  to 
0(M  log  M).  However,  we  consider  this  as  a  preliminary  result  since  it  is  also 
important  to  assure  that  the  break  even  point  with  the  usual  dense  matrix- 
vector  multiply  is  low  and  the  fast  scheme  becomes  beneficial  for  a  relatively 
small  M.  We  note  that  the  overall  operation  count  also  takes  into  account  the 
number  of  iterations  Nlt  for  solving  the  nonlinear  system  (11)  in  Section~3.3. 
Since  the  number  of  iterations  depends  on  the  length  of  the  time  interval  on 
which  we  discretize  our  system  of  ODEs,  it  may  grow  as  a  function  of  M  and 
diminish  the  benefit  of  increasing  the  number  of  nodes  M. 

We  now  outline  two  approaches  for  applying  the  integration  matrix  S  in 
O(MlogM)  operations. 

5.1  The  Partitioned  Low  Rank  (PLR)  representation 

The  Partitioned  Low  Rank  (PLR)  representation  may  be  applied  in  a  variety 
of  problems  where  the  off-diagonal  part  of  a  matrix  has  a  relatively  small  rank. 
In  particular,  PLR  may  be  used  for  integration  and  differentiation  matrices. 
The  idea  is  to  subdivide  the  matrix  as  in  Figure  3  and  then  represent  individual 
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off-diagonal  blocks  as  a  sum  of  rank  one  matrices, 


J2 

k= 1 

where  and  v*.,  k  —  1, . . . ,  r  are  vectors  of  appropriate  size  for  a  given  block. 
In  this  representation,  the  number  of  terms  r  (the  rank  of  the  off-diagonal 
block)  is  selected  for  a  given  user-supplied  accuracy  and  may  be  found  by 
the  Singular  Value  Decomposition  (SVD).  However,  since  we  do  not  require 
orthogonality  between  vectors,  a  simpler  algorithm  may  be  used  instead.  If 
the  ranks  of  off-diagonal  blocks  r  are  small,  the  cost  of  applying  matrces  in 
the  PLR  representation  is  0(M  log  M).  For  more  details,  see  [2,5]. 


Figure  3.  Partitioning  of  a  matrix  in  PLR  representation 

As  an  example,  consider  the  integration  matrix  in  Theorem  4  with  70  nodes, 
a  size  that  we  have  found  to  give  a  good  compromise  between  large  number  of 
nodes  vs.  small  number  of  iterations.  In  this  case,  optimal  PLR  performance 
is  obtained  by  using  just  one  level  of  subdivision  with  the  rank  r  =  13  for  the 
lower  and  r  =  14  for  the  upper  off-diagonal  blocks  to  attain  double  precision 
accuracy.  In  this  case,  applying  the  integration  matrix  in  the  PLR  representa¬ 
tion  requires  7735  operations,  whereas  applying  the  dense  integration  matrix 
takes  9800  operations,  only  a  moderate  improvement  in  performance. 

It  remains  a  challenge  to  increase  the  intervals  of  integration  (in  order  to  work 
with  a  larger  matrix  size  M )  without  increasing  the  number  of  iterations. 

5.2  Quadratures  of  Gauss- Trapezoidal  type 

hi  [1]  Alpert  introduced  high  order,  polynomial  based  quadratures  with  nodes 
equally  spaced  in  the  interior  of  the  interval  and  with  a  limited  region  of 
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higher  node  density  near  its  boundary.  Such  node  distribution  allows  one  to 
use  Fourier-type  methods  combined  with  a  low  rank  correction.  It  turns  out 
that  a  similar  quadrature  for  exponentials  (bandlimited  functions)  may  also 
be  constructed  and  we  use  an  example  of  such  quadratures  (supplied  to  us 
by  Brad  Alpert)  in  our  numerical  experiments.  We  note  that  this  arrange¬ 
ment  of  nodes  does  reduce  the  bandlimit  of  functions  accurately  integrated 
by  the  quadrature  but  its  structure  provides  certain  advantages  vis-a-vie  fast 
algorithms. 

Let  us  demonstrate  how  this  quadrature  may  be  used  for  fast  application  of 
the  integration  matrix  defined  by 

rSk 

Ski  =  J  |  Ri(s)  ds. 

We  note  that  although  currently  we  do  not  have  a  proof  that  this  matrix 
is  symplectic,  we  verified  numerically  that  it  is  very  close  to  the  symplectic 
integration  matrix  defined  in  Theorem  4. 


Using  representation  of  the  interpolating  functions  Ri(s)  via  exponentials  (i.e. , 
their  definition  in  |5]),  we  have 


M 


icOjS 


and  write  Ski  as 
where 

and 


Ri(s )  =  E  AUe 

3= 1 


Ski  —  Ski  —  oil 

M  eic9kGj 

Ski  =  J2Aj— 


(16) 


at  =  Y^A 


j=l  ic6J 


M  g— icdj 


iy 


3  = 1 


icdj 


We  note  that  matrix  elements  cq  are  independent  of  the  index  k  and,  therefore, 
may  be  applied  as  a  rank-one  matrix.  Hence,  we  only  need  to  consider  applying 
S.  Let  us  define 

Ekj  =  eic6kdl 

and  the  diagonal  matrix 


Dkj 


ik’i  =  j 

0  i  ^  j 


Since  functions  Ri(s)  are  interpolating,  Ri(6k)  =  Ski,  and  we  have  from  (16) 
that  A  =  E~l .  Hence,  we  write 


S  =  E~lDE 
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Let  us  introduce  matrix  of  the  disrete  Fourier  Transform, 


TP,  _  2-nikj/M 

A  KJ  °  1 

and  a  parameter 

L2ti 
C  =  2 M’ 

where  L  =  2a  +  ninterior  ~  1,  '^interior  is  the  number  of  interior  equispaced  nodes 
and  a  is  the  size  of  the  region  near  the  boundary  with  nodes  of  higher  density. 
We  then  have  for  the  interior  nodes 

—  eicdk6i  _  eic(2k/L-\)(2j/L-l)  _  ^2mkj/M^—2ick/L^—2icj/L^ic 

or 

Ekj  =  Fkje-2ick/Le~2icj/Leic, 

where 

Df  =  diag  {e-^'/V/2}^1 

For  the  full  matrix  E  we  write 

E  =  DfFDf  +  G 

where  G  is  a  matrix  that  with  zeros  in  its  interior  and  a  non-zero  border  near 
the  edges  of  the  matrix.  Thus,  matrix  G  has  low-rank. 

Applying  matrix  S'  to  a  vector  /  may  now  be  written  as 

ST  =  ( DfFDf  +  G)-1  (D{DFFDF  +  G)f)  =  ( DfFDf  +  G)-1g, 

where 

g  =  ( D(DfFDf  +  G)f) 

We  note  F  is  a  matrix  of  the  Discrete  Fourier  Transform  (DFT)  which  may 
be  applied  in  0(M  log  M)  operations  using  the  FFT.  Since  G  is  of  low  rank 
and  D  and  DF  are  diagonal,  g  may  be  computed  in  0(M  log  M)  operations. 
Furthermore,  we  have 

(. DfFDf  +  G)~l  =  (/  +  Df~1F~1Df~1G)~1F~1 

and,  using  the  fact  that  Df,^1F~1Df~1G  is  of  low  rank,  we  have  a  fast  algo¬ 
rithm  to  compute  Sf  by  the  use  low  rank  update  of  the  inverse. 
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6  Numerical  Examples 


6.1  Orbit  Determination 

Current  ODE  integrators  for  orbit  computations  use  high  order  Gauss-Jackson 
multistep  methods.  Several  years  ago  we  have  compared  the  number  of  func¬ 
tion  evaluations  of  this  method  with  an  explicit  version  of  the  deferred  spec¬ 
tral  correction  method  in  [7,9].  The  comparison  was  in  favor  of  the  latter 
and  it  showed  that  high  order  multistep  methods  are  significantly  oversam¬ 
pled  (which  also  follows  from  theoretical  considerations) .  Therefore,  a  possible 
gain  in  performance  should  come  from  the  reduced  sampling  requirements;  we 
believe  that  quadratures  for  bandlimited  functions  are  the  right  tool  for  this 
purpose.  The  new  methods  will  add  to  the  variety  of  available  techniques  and 
will  also  provide  new  time  evolution  schemes  for  partial  differential  equations. 

Let  us  consider  the  spherical  harmonic  model  of  the  gravitational  potential, 
V(x,y,z),  which  in  the  spherical  system  of  coordinates  is  written  in  terms  of 
the  spherical  harmonics, 


(17) 


Let  G  =  (Gx,  Gy ,  GZY  =  W  and  we  use  the  gravity  model  of  the  form 


(x(t)2+y(t)2+z{t)2)3/2 


In  order  to  reduce  the  number  of  function  calls  to  the  full  model,  we  initially 
use  the  gravity  model  with  N  —  3  on  a  large  portion  (e.g.,  1/2  of  a  period)  of 
an  orbit  to  solve  the  system  of  nonlinear  equations  via  a  fixed  point  iteration. 
We  then  access  the  full  gravity  model  and  evaluate  the  gravitational  force  at 
the  nodes  that  by  now  are  located  close  to  their  correct  positions.  We  continue 
iteration  (without  accessing  the  full  gravity  model  again)  to  adjust  the  orbit. 
This  results  in  an  essentially  correct  trajectory.  At  this  point  we  may  (and 
currently  do)  access  the  full  gravity  model  one  more  time  to  evaluate  the 
gravitational  force  and  perform  another  iteration.  Thus,  we  access  the  full 
gravity  model  at  most  twice  per  node  while  using  a  number  of  nodes  that  is 
substantially  lower  than  in  traditional  methods. 

Let  us  now  re-write  the  orbit  determination  problem  (using  only  gravitational 
forces)  in  a  form  that  conforms  with  the  algorithm  in  Section  3.3.  We  define 
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the  six  component  vector 


y  (t)  = 


x(t) 

x'(t) 

y(t) 

y'(t ) 
z(t) 
z'[t) 


where  [x,  y,z]*  denote  the  positions,  and  \x'  ,y\z'Y  denote  the  velocities.  We 
also  define  the  matrix 


010000 

000000 

000100 

000000 

000001 

000000 


L  = 


and  the  right  hand  side 


g(f>y(f))  = 


'  (x(t)2+y(t)2+z(t)2)3/2  +  En=2  Gx(y(t)) 
0 

'  (x{t)2+y(t)2+z(t)2)3/2  +  Sn=2  Gy(y(t)) 


(x(t) 


2+y(^)2|2(t)2)3/2  +  T,n= 2  G*(y(t))  _ 


where  G(y(£))  =  [Gx(y(t)),  Gy(y(t)),  Gz(y(i))]  denotes  the  x,  y  and  ^-components 
of  gravity  model  starting  with  n  =  2  and  up  to  order  N .  The  ODE  describing 
the  orbit  determination  problem  is  now  given  by  (10). 

We  apply  the  algorithm  in  Section  3.3  by  first  using  Nlt  =  2  and  N  =  3 
to  obtain  an  approximate  solution.  We  then  switch  to  the  full  model  with 
IV  =  70  and,  at  this  point  require  only  one  or  two  iterations  Nit.  Thus,  we 
need  to  access  the  full  model  only  once  or  twice  per  node. 
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6.2  Timing 


We  tested  and  timed  the  algorithm  by  computing  an  orbit  of  86000  seconds 
(~  1  day)  using  a  36  and  70  degree  gravitational  model.  The  timing  was 
performed  on  a  computer  with  an  Intel  Core2  Extreme  processor  at  2.96  GHz 
and  4  GB  of  RAM  at  800  MHz.  (However,  only  a  small  fraction  of  the  RAM 
was  actually  used  during  the  computation.) 

6.2.1  Timing  results  using  spherical  harmonics 

We  Erst  timed  the  algorithm  using  the  spherical  harmonic  gravitational  model. 
The  timing  results  are  given  in  Table  1. 


Degree 

CPU  Time  (s) 

Number  of  function  calls 

36 

9.0E-02 

6160 

70 

2.7E-01 

6160 

Table  1 

CPU  time  and  number  of  function  calls  to  full  gravitational  model  when  computing 
a  orbit  of  86000  seconds  using  spherical  harmonics. 


6.2.2  Timing  results  using  a  cubed  sphere  spline  model 

We  also  timed  the  algorithm  when  computing  the  gravitational  model  using 
a  local  cubed  sphere  spline  model  [3].  The  timing  results  are  given  in  Table  2. 


Degree 

CPU  Time  (s) 

Number  of  function  calls 

41 

3.2E-02 

6160 

70 

3.2E-02 

6160 

Table  2 

CPU  time  and  number  of  function  calls  to  full  gravitational  model  when  computing 
a  orbit  of  86000  seconds  using  a  local  cubed  sphere  model. 


6.2.3  Current  (indirect,  crude)  estimate  of  timing 

We  tested  the  speed  of  our  code  on  a  sample  NOAA  satellite  orbit,  for  which 
we  received  timing  results  using  the  SPEPH  code.  Since  the  SPEPH  code  and 
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our  code  ran  on  different  machines,  we  were  able  only  to  estimate  the  timing 
using  the  following  approach.  We  were  given  the  speed  of  SPEPH  code  for 
gravitational  models  with  N  =  36  and  N  =  70  order  and  degree.  Assuming 
that  timing  captures  the  cost  of  integration  and  access  to  the  gravity  model 
and  assuming  that  the  cost  of  integration  is  the  same  for  both  models,  we 
write 

Adeg)  _  ,  .(deg) 

Ltotal  —  L code  “r  l modeli 

where  tcode  is  the  time  of  integration  and  t^Jel  is  the  time  of  access  to  the 
model.  We  further  assume  that 

=  PO/36)2  ss  3.78, 

since  the  cost  of  access  to  the  spherical  harmonic  model  grows  quadratically 
with  its  order  and  degree.  Given  this  information  for  SPEPH  code  and  timing 
our  code  we  obtain  the  following  comparison  by  solving  the  above  equations 


t 


(36) 

model 


,(70) 

k  total 


t 


(36) 

total 


3.78  -  1 


n  36  •  (tm  -  t(36)  1 

U.OU  \itotai  i  total) 


,(70)  _  q  70  .  ^ 

"  model  O.IO  Lmol}el 


+  _  ,(36) 

Lcode  ‘'total 


t. 


(36) 
modt 
(36) 
model  * 


1  36  •  —  P'30-’  ) 

l.OU  \^total  ^ total ) 


t 

(36) 


Using  the  values  from  timing  SPEPH,  tf?al  =  5.21  and  t(7Jl  =  16.36,  we  get 


,(36)  _ 

Lmodel  ^ 


^ model 
t code  ^ 


t 


(70) 

model 

tcode 


For  our  code  we  have  tf^al  =  6.3  and  t\'^’al  =  18.6,  and  we  get 


4.0 

15.12 

1.21 

a  12.5 

(70) 


(36) 

Lmodel 

470) 


4.4 
16.63 

tcode  ~  1-9 


tK  1 

model 


t 


(70) 

model 


tcode 


8.75 


Thus,  currently,  our  implicit  code  appears  to  be  about  30%  slower  than  the 
explicit  SPEPH  code  (if  the  indirect  estimate  above  is  correct).  Also,  it  appears 
that  the  computation  of  the  spherical  harmonics  has  been  accelerated  in  the 
current  SPEPH  code  in  comparison  with  its  predecessor  since  the  tests  of  6  — 7 
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years  ago  were  done.  For  this  reason,  the  local  gravity  model  developed  in  [3] 
is  faster  than  the  spherical  harmonic  model  by  a  lower  factor,  perhaps  only  a 
factor  of  10  rather  than  40  as  before.  However,  this  may  be  changed  to  higher 
performance  factors  so  that  the  local  model  becomes  faster  since  it  is  a  simple 
matter  of  trading  speed  vs  memory  (memory  is  no  longer  a  problem  for  the 
sizes  needed  for  these  models).  We  also  need  to  accelerate  our  code  to  match 
and,  hopefully,  exceed  the  speed  of  the  current  SPEPH  code.  Note  that  the 
our  code  is  A-stable  and  symplectic  (if  need  be).  The  code  is  currently  being 
tested  for  accuracy  by  Brandon  Jones  and  appears  to  match  the  results  of  the 
codes  used  in  CU  Aerospace  Department. 


7  Conclusions 

We  have  constructed  an  implicit,  symplectic  integrator  that  has  speed  compa¬ 
rable  to  the  explicit  multistep  integrator  currently  used  for  orbit  computation. 
In  combination  with  a  local  model  for  gravity  it  achieves  a  factor  of  8.3  im¬ 
provement  in  speed  compared  to  using  a  global  gravity  model. 
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